x_mesh:
[[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]
[0. 1. 2. 3. 4. 5. 6. 7. 8. 9.]]
y_mesh:
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[2. 2. 2. 2. 2. 2. 2. 2. 2. 2.]
[3. 3. 3. 3. 3. 3. 3. 3. 3. 3.]
[4. 4. 4. 4. 4. 4. 4. 4. 4. 4.]
[5. 5. 5. 5. 5. 5. 5. 5. 5. 5.]
[6. 6. 6. 6. 6. 6. 6. 6. 6. 6.]
[7. 7. 7. 7. 7. 7. 7. 7. 7. 7.]
[8. 8. 8. 8. 8. 8. 8. 8. 8. 8.]
[9. 9. 9. 9. 9. 9. 9. 9. 9. 9.]]
Tutorial 9: Advanced Modelling and Visualisation
Programming for Engineers (ENGG1001)
Tutorial Tasks
Task 1: Creating a Meshgrid for 3D Surface Plots
So far, you have had practice with creating 1D numpy arrays and plotting them. For 3D surface plots, we need to create a grid of points in the x-y plane. Write a function create_meshgrid(length: float, num_points: int) -> tuple[np.ndarray, np.ndarray] that creates a meshgrid of points from -length to length with num_points for both x and y axes. For a further explanation of what a meshgrid is, see the callout below.
What exactly is a meshgrid? A meshgrid is a 2D grid of points that can be used for plotting surfaces or evaluating functions over a 2D domain. It consists of two 2D arrays. Consider the example below where we create a 10 x 10 meshgrid:
Essentailly, a meshgrid is a way of presenting a set of cartesian coordinates in a structured format. So if we feed an index (which essentially acts as a coordinate) into the meshgrid, we can get the corresponding x and y values. For example, if we want to get the x and y values at the index (2, 3), we can do:
x value at x_mesh[2, 3]: 3.0
y value at y_mesh[2, 3]: 2.0
Here is a sample usage of the function:
What function could you use to create the 1D arrays for x and y? Once you have created the 1D arrays, then you can use numpy.meshgrid to create the 2D grid.
meshgrid.py
import numpy as np
def create_meshgrid(length: float, num_points: int) -> np.ndarray:
"""Creates a meshgrid of points from -length to length with num_points for both x and y axes.
Parameters
----------
length : float
The half-length of the grid in both x and y directions.
num_points : int
The number of points to generate along each axis.
"""
x = np.linspace(-length, length, num_points)
y = np.linspace(-length, length, num_points)
x_mesh, y_mesh = np.meshgrid(x, y)
return x_mesh, y_mesh
x_mesh, y_mesh = create_meshgrid(10, 100)
print(x_mesh.shape)
print(x_mesh[0, 0])
print(y_mesh[50, 0])(100, 100)
-10.0
0.10101010101010033
Task 2: Creating a Plane Surface
Now that we can create a meshgrid of points, we can use it to create a 3D surface. Write a function plane_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, a: float, b: float, c: float) -> np.ndarray that computes the z values for a plane surface defined by the equation: \[z = ax + by + c\]
Here is a sample usage of the function:
To compute the wave surface, you will need to calculate the radial distance R for each point in the meshgrid and then apply the RBF formula. Remember to use vectorisation!
plane_surface.py
def plane_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, a: float, b: float, c: float) -> np.ndarray:
"""Computes the z values for a plane surface defined by the equation z = ax + by + c.
Parameters
----------
x_mesh : np.ndarray
The meshgrid of x coordinates.
y_mesh : np.ndarray
The meshgrid of y coordinates.
a : float
The coefficient for x.
b : float
The coefficient for y.
c : float
The constant term.
Returns
-------
np.ndarray
The computed z values for the plane surface.
"""
z = a * x_mesh + b * y_mesh + c
return z
x_mesh, y_mesh = create_meshgrid(10, 100)
plane = plane_surface(x_mesh, y_mesh, 1.0, 2.0, 3.0)
print(plane.shape)
print(plane[0, 0])
print(plane[50, 50])(100, 100)
-27.0
3.303030303030301
Task 3: Modelling waves with a Radial Basis Function
Now that we have a meshgrid of points, we can model a wave surface using a Radial Basis Function (RBF). The RBF is defined as: \[f(x, y) = A \cdot \cos(kr - \omega t) \cdot e^{-\alpha r}\] Where:
- \(A\) is the amplitude of the wave
- \(k\) is the wave number
- \(r\) is the radial distance from the origin: \(r = \sqrt{x^2 + y^2}\)
- \(\omega\) is the angular frequency
- \(\alpha\) is the decay constant of the wave with respect to distance
- \(t\) is time in seconds
Write a function radial_wave(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float, k: float, omega: float, alpha: float, t: float) -> np.ndarray that computes the wave surface values for the given parameters. Here is a sample usage of the function:
radial_wave.py
def radial_wave(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float, k: float, omega: float, alpha: float, t: float) -> np.ndarray:
"""Computes the wave surface values for the given parameters using a Radial Basis Function.
Parameters
----------
x_mesh : np.ndarray
The meshgrid of x coordinates.
y_mesh : np.ndarray
The meshgrid of y coordinates.
A : float
The amplitude of the wave.
k : float
The wave number.
omega : float
The angular frequency.
alpha : float
The decay constant of the wave with respect to distance.
t : float
Time in seconds.
Returns
-------
np.ndarray
The computed wave surface values.
"""
r = np.sqrt(x_mesh**2 + y_mesh**2)
wave_surface = A * np.cos(k * r - omega * t) * np.exp(-alpha * r)
return wave_surface
x_mesh, y_mesh = create_meshgrid(10, 100)
wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)
print(wave_surface.shape)
print(wave_surface[50, 50])
print(wave_surface[0, 0])(100, 100)
0.9458561801297212
-0.24310473049518086
Task 4: Visualising 3D Surfaces
Now that we have a function to compute the wave surface, we can visualise it using Matplotlib’s 3D plotting capabilities. Write a function plot_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, surface: np.ndarray) that creates a 3D surface plot of the surface.
Here is a sample usage of the function potting the plane surface:
Which produces the following 3D surface plot:
Also try plotting the wave surface using the same function:
Which produces the following 3D surface plot:
To be able to setup an axis for a 3D plot, you need to specify the projection argument when creating the subplot. The plot_surface method can be used to create the surface plot. You can use the code below to get started:
If you want to make your plot look nicer, you can also specify a colormap using the cmap argument in plot_surface. For example, you can use cmap='viridis' to use the Viridis colormap, like the example above.
plot_surface.py
import matplotlib.pyplot as plt
def plot_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, surface: np.ndarray):
"""Creates a 3D surface plot of the surface.
Parameters
----------
x_mesh : np.ndarray
The meshgrid of x coordinates.
y_mesh : np.ndarray
The meshgrid of y coordinates.
surface : np.ndarray
The computed surface values.
"""
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
ax.plot_surface(x_mesh, y_mesh, surface, cmap='viridis')
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Surface')
ax.set_title('3D Surface Plot')
plt.show()
x_mesh, y_mesh = create_meshgrid(10, 100)
wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)
plot_surface(x_mesh, y_mesh, wave_surface)Now that you have created an easy way to visualise the wave surface, determine what effects the following parameters have on the wave surface by changing their values and observing the resulting plot:
- Amplitude
A - Wave number
k - Angular frequency
omega - Decay constant
alpha - Num points in the meshgrid
num_points
Task 5: identifying local maxima and minima
Now that we have a visualisation of the wave surface, we can identify the local maxima and minima of the surface. Write a function find_local_extrema(surface: np.ndarray) -> tuple[np.ndarray, np.ndarray] that finds and returns the radius values (how far from the origin) of the local maxima and minima of the surface. Here is a sample usage of the function:
To find the local maxima and minima of radial function, we can use the fact that the surface is “radially symmetric”, meaning that the surface values only depend on the radius from the origin. As such, we can compute the radius values for a central “slice” of the surface (e.g. a given row or column of the surface), and then find the local maxima and minima of that slice.
To illustrate this, consider the following central slice of the curve above:
We can see that the local maxima and minima of the surface correspond to the local maxima and minima of this slice. Therefore, we can find the local maxima and minima of this slice to determine the radius values of the local maxima and minima of the surface.
To find the local maxima and minima of a 1D array, you can start by using the np.diff function to compute the differences between consecutive elements. Then you can find where the sign of these differences changes (i.e. from positive to negative for maxima, and from negative to positive for minima). Then you can store the corresponding radius values for those indices and return them as the output of the function.
Two solutions are presented below, one using a more manual approach, and the other using vectorisation:
find_local_extrema.py
def find_local_extrema_simple(surface: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Finds and returns the radius values of the local maxima and minima of the surface.
Parameters
----------
surface : np.ndarray
The computed surface values.
Returns
-------
tuple[np.ndarray, np.ndarray]
A tuple containing two arrays: the first array contains the radius values of the local maxima, and the second array contains the radius values of the local minima.
"""
# Compute the radius values for a slice of the
# surface (e.g. the middle row)
radius = np.sqrt(x_mesh[x_mesh.shape[0] // 2]**2
+ y_mesh[y_mesh.shape[0] // 2]**2)
diffs = np.diff(surface[surface.shape[0] // 2]) # taking the middle row as the slice
maxima = []
minima = []
for i in range(1, len(diffs) - 1):
if diffs[i-1] > 0 and diffs[i] <= 0: # local maximum
maxima.append(radius[i])
elif diffs[i-1] < 0 and diffs[i] >= 0: # local minimum
minima.append(radius[i])
return np.array(maxima), np.array(minima)
def find_local_extrema(surface: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
"""Finds and returns the radius values of the local maxima and minima of the surface.
Parameters
----------
surface : np.ndarray
The computed surface values.
Returns
-------
tuple[np.ndarray, np.ndarray]
A tuple containing two arrays: the first array contains the radius values of the local maxima, and the second array contains the radius values of the local minima.
"""
# Compute the radius values for a slice of the
# surface (e.g. the middle row)
radius = np.sqrt(x_mesh[x_mesh.shape[0] // 2]**2
+ y_mesh[y_mesh.shape[0] // 2]**2)
# Compute the differences between consecutive elements in the surface slice
diff = np.diff(surface[surface.shape[0] // 2]) # taking the middle row as the slice
# Find where the sign of the differences changes to identify local maxima and minima
maxima_indices = np.where((diff[:-1] > 0) & (diff[1:] <= 0))[0] + 1
minima_indices = np.where((diff[:-1] < 0) & (diff[1:] >= 0))[0] + 1
# Get the corresponding radius values for those indices
maxima = radius[maxima_indices]
minima = radius[minima_indices]
return maxima, minima
x_mesh, y_mesh = create_meshgrid(10, 100)
wave_surface = radial_wave(x_mesh, y_mesh, 1.0, 2.0, 1.0, 0.1, 0.0)
maxima, minima = find_local_extrema(wave_surface)
print(maxima)
print(minima)[9.39448244 6.16244406 3.1329419 0.14284985 3.1329419 6.16244406
9.39448244]
[7.77843366 4.7485492 1.51851479 1.51851479 4.7485492 7.77843366]
Task 6: Animating a 3D Wave Surface
For a challenge task, try to create an animation of the wave surface over time. However there is a slight twist: the equation for the wave surface now accounts for time decay as well, so the wave surface will not only change with time, but it will also decay over time. The new equation is: \[f(x, y, t) = A \cdot \cos(kr - \omega t) \cdot e^{-\alpha r} \cdot e^{-\beta t}\] Where \(\beta\) is the decay constant of the wave with respect to time. You can use Matplotlib’s animation capabilities to create the animation. Here is a sample usage of the function:
Use the following incomplete code as a starting point for your animation function. You will need to complete the update_image function and the animate_wave_surface function to create the animation.
from matplotlib import animation
def update_image(frame, x_mesh, y_mesh, A, k, omega, alpha, beta, ax):
"""Updates the wave surface for each frame of the animation.
Parameters
----------
frame : int
The current frame number.
x_mesh : np.ndarray
The meshgrid of x coordinates.
y_mesh : np.ndarray
The meshgrid of y coordinates.
A : float
The amplitude of the wave.
k : float
The wave number.
omega : float
The angular frequency.
alpha : float
The decay constant of the wave with respect to distance.
beta : float
The decay constant of the wave with respect to time.
ax : matplotlib.axes._subplots.Axes3DSubplot
The 3D axis object to update.
"""
wave_surface = # write your code here!
ax.clear()
ax.plot_surface(x_mesh, y_mesh, wave_surface)
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Wave Surface')
ax.set_title('Animated 3D Wave Surface Plot')
return [ax]
def animate_wave_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float,
k: float, omega: float, alpha: float, beta: float,
delta_t: float, num_steps: int):
"""Creates an animation of the wave surface over time.
Parameters
----------
x_mesh : np.ndarray
The meshgrid of x coordinates.
y_mesh : np.ndarray
The meshgrid of y coordinates.
A : float
The amplitude of the wave.
k : float
The wave number.
omega : float
The angular frequency.
alpha : float
The decay constant of the wave with respect to distance.
beta : float
The decay constant of the wave with respect to time.
delta_t : float
The time step for each frame of the animation.
num_steps : int
The number of steps in the animation.
"""
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
#| decorate the plot with labels and title here!
anim = animation.FuncAnimation(fig, update_image, frames=np.arange(0, num_steps, delta_t),
fargs=(x_mesh, y_mesh, A, k, omega, alpha, beta, ax), interval=50)
plt.show()animate_wave_surface.py
from matplotlib import animation
def update_image(frame, x_mesh, y_mesh, A, k, omega, alpha, beta, ax):
"""Updates the wave surface for each frame of the animation.
Parameters
----------
frame : int
The current frame number.
x_mesh : np.ndarray
The meshgrid of x coordinates.
y_mesh : np.ndarray
The meshgrid of y coordinates.
A : float
The amplitude of the wave.
k : float
The wave number.
omega : float
The angular frequency.
alpha : float
The decay constant of the wave with respect to distance.
beta : float
The decay constant of the wave with respect to time.
ax : matplotlib.axes._subplots.Axes3DSubplot
The 3D axis object to update.
"""
R = np.sqrt(x_mesh**2 + y_mesh**2)
wave_surface = A * np.cos(k * R - omega * frame) * np.exp(-alpha * R) * np.exp(-beta * frame)
ax.clear()
ax.plot_surface(x_mesh, y_mesh, wave_surface)
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Wave Surface')
ax.set_title('Animated 3D Wave Surface Plot')
return [ax]
def animate_wave_surface(x_mesh: np.ndarray, y_mesh: np.ndarray, A: float,
k: float, omega: float, alpha: float, beta: float,
delta_t: float, num_steps: int):
"""Creates an animation of the wave surface over time.
Parameters
----------
x_mesh : np.ndarray
The meshgrid of x coordinates.
y_mesh : np.ndarray
The meshgrid of y coordinates.
A : float
The amplitude of the wave.
k : float
The wave number.
omega : float
The angular frequency.
alpha : float
The decay constant of the wave with respect to distance.
beta : float
The decay constant of the wave with respect to time.
delta_t : float
The time step for each frame of the animation.
num_steps : int
The number of steps in the animation.
"""
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
ax.set_xlabel('X-axis')
ax.set_ylabel('Y-axis')
ax.set_zlabel('Wave Surface')
ax.set_title('Animated 3D Wave Surface Plot')
anim = animation.FuncAnimation(fig, update_image, frames=np.arange(0, num_steps, delta_t),
fargs=(x_mesh, y_mesh, A, k, omega, alpha, beta, ax), interval=50)
plt.show()